// Gamma function
// Copyright 2020 Glenn McIntosh
// licensed under the GNU General Public Licence version 3
// uses Lanczos approximation
// N is number of terms in series, must be even
#pragma once
#include <cmath>
#include <complex>
template<int N>
class Gamma
{
public:
	Gamma()
	{
		// odd chebyshev polynomials
		int64_t cheby[N+1][N+1] {};
		cheby[0][0] = 1;
		for (int n = 1; n <= N; ++n)
			cheby[n][0] = cheby[n][0]*2 - cheby[n-1][0];
		int64_t b1 = 1;
		for (int l = 1; l <= N; ++l)
		{
			for (int n = l; n <= N; ++n)
			{
				b1 = cheby[n-1][l-1]*2 - b1;
				cheby[n][l] = b1*2 - cheby[n-1][l];
			}
			b1 = 0;
		}

		// multiply by binomial and sloan series
		for (int n = 1; n <= N; ++n)
			for (int l = 0; l <= N; ++l)
				cheby[0][l] += cheby[n][l];
		int64_t sloan = -1;
		for (int n = 1; n <= N; ++n)
		{
			int64_t b{-(n+n)};
			for (int j = n+1; j <= N; ++j)
			{
				for (int l = 0; l <= N; ++l)
					cheby[n][l] += b * cheby[j][l];
				b *= -(j+n); b /= j-n+1;
			}
			for (int l = 0; l <= N; ++l)
				cheby[n][l] *= sloan;
			sloan *= (2*n)*(2*n+1);
			sloan /= n*n;
		}

		// prescale
		int64_t ff{1};
		for (int l = 0; l <= N; ++l)
		{
			for (int n = 0; n <= N; ++n)
			{
				cheby[n][l] /= 1<<l;
				cheby[n][l] *= ff*2;
			}
			ff *= l*2+1;
		}
		cheby[0][0] /= 2;

		// multiply factors
		for (int n = 0; n <= N; ++n)
		{
			long double sum{};
			for (int l = 0; l <= N; ++l)
			{
				long double t = exp(g2 + l) / sqrt(pow(g2+l, l*2+1));
				sum += cheby[n][l] * t;
			}
			coeff[n] = sum;
		}
	}
	template<typename T> T operator()(T s)
	{
		bool invert = std::real(s) < 0.5;
		if (invert)
			s = 1.-s;
		T y = s;
		T sum = coeff[0];
		for (int n = 1; n <= N; ++n)
			sum += coeff[n] / (y + n*1.);
		sum /= s;
		sum *= exp(s*log(s+g2*1.)) * sqrt(s+g2*1.) / exp(s+g2*1.);
		if (invert)
			sum = tau/2 / (sin(tau/2*s) * sum);
		return sum;
	}
private:
	static constexpr int g2{6}; // g+0.5
	static constexpr double tau = ::acos(-1)*2;
	double coeff[N+1] {};
};
